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ABSTRACT 

Astrophysical relativistic flow problems require high resolution three-dimensional numerical simu- 
lations. In this paper, we describe a new parallel three-dimensional code for simulations of special 
relativistic hydrodynamics (SRHD) using both spatially and temporally structured adaptive mesh 
refinement (AMR). We used the method of lines to discretize the SRHD equations spatially and 
a total variation diminishing (TVD) Runge-Kutta scheme for time integration. For spatial recon- 
struction, we have implemented piecewise linear method (PLM), piecewise parabolic method (PPM), 
third order convex essentially non-oscillatory (CENO) and third and fifth order weighted essentially 
non-oscillatory (WENO) schemes. Flux is computed using either direct flux reconstruction or approx- 
imate Riemann solvers including HLL, modified Marquina flux, local Lax-Friedrichs flux formulas and 
HLLC. The AMR part of the code is built on top of the cosmological Eulerian AMR code enzo. We 
discuss the coupling of the AMR framework with the relativistic solvers. Via various test problems, 
we emphasize the importance of resolution studies in relativistic flow simulations because extremely 
high resolution is required especially when shear flows are present in the problem. We also present 
the results of two 3d simulations of astrophysical jets: AGN jets and GRB jets. Resolution study of 
those two cases further highlights the need of high resolutions to calculate accurately relativistic flow 
problems. 

Subject headings: hydrodynamics-methods: numerical-relativity 



1. INTRODUCTION 

Relativistic flow problems are important in many 
astrophysical phenomena including gamma-ray burst 
(GRB), active galactic nuclei (AGN), as well as micro- 
quasar and pulsar wind nebulae, among others. Appar- 
ent super luminal motion is observed in many jets of ex- 
tragalactic radio sources associated with AGN. Accord- 
ing to the currently accepted standard model, this im- 
plies the jet flow velocities as large as 99% of the spee d 
of light (jBlandford et al.lll977l : [Begelman et al. Ill984f ). 
Similar phenome na are also seen in microquas ars such as 
GRS 191 5+105 (iMirabel fc Rodriguez I fl99l and GRO 
J 1655-40 ([Tingay et al. Ifl995l ) thus by similar arguments 
relativistic flows are thought to play a role. In the 
case of GRB, the observed non-thermal spectrum im- 
plies that the source must be optically thin, which can 
be used to put a lim it on the minimum L orentz factor 
within those bursts (|Lithwick fc Saril l200lh . This argu- 
ment shows that the source of GRB must be in highly rel- 
ativistic motion. This conclusion is further confirmed by 
the rapidly increasing number of GRB afterglow o bser- 
vations by the Swift satellite (Ge hrels et al. 1 12004). To 
understand physical processes in those phenomena quan- 
titatively, high resolution multi-dimensional simulations 
are crucial. 

Jim Wilson and collaborators pioneered the numerical 



solut i on of relativistic hydrodynamics equa tions ( Wilson! 
119721 : ICentrella fc Wilsonl 119841 : lHawlev et all Il984al lrl). 
Starting with these earliest papers this has typically been 
done in the context of general relativistic problems such 
as accretion onto black holes and super novae explosions. 
The problem was recognized to be diffi cult to solve when 
the L orentz factor becomes large (|Norman fc Wi nkler 
Il986f ) and a solution with an implicit adaptive scheme 
was demonstrated in one dimension. Unfortunately, this 
approach is not generalizable to multi-dimensions. How- 
ever, in the past two decades accurate solvers based on 
Godunov's scheme have been designed that have adopted 
shock-capturing schemes for Newtonian fluid to the rela- 
tivist ic fluid equations in conservation form (for a review 
see iMartf fc Mulleri [2003). Such schemes, called high- 
resolution shock-capturing (HRSC) methods, have been 
proven to be very useful in capturing strong discontinu- 
ities with a few numerical zones without serious numer- 
ical oscillations. We will discuss a number of them in 
section 3. 

Studies involving astrophysical fluid dynamics in gen- 
eral are benefiting tremendously from using spatial 
and temporal adapti ve techniques. Smoothed parti- 
cle h ydrodynamics (|Gingold fc Monaghanl Il977l : iLucyl 
Il977f ) is a classic example by bei ng; a Lagrangian 
meth od. Increasingly also variants of iBerger fc C olella 



([I989[ )'s structured adaptive mesh technique are be- 
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tic hydrodynamics (iHughes et al. I2QQ2I; lAnninos et al.l 



2005; Zhang & MacFadvenl 120061 : iMorsony et al. 1 12006: 
Meliani et al. I 120071 ) where certainly the work of 
Hughes et al.l (l2QQ2h showed that a serial AMR code 



could solve problems even highly efficient parallel fixed 
grid codes would have difficulty with. In this paper we 
discuss our implementation of different hydrodynamics 
solvers with various reconstruction schemes as well as 
different time integrators on top of th e enzo framework 
previo usl y developed for co smology ( Bryan fc Normanl 
Il997al lbl: [Bryan et all l200ll : lO'Shea et alJl2QQ4f i. This 
new code we call renzo is adaptive in time and space 
and is a dynamically load balanced parallel using the 
standard message passing interface. 

In the following we briefly summarize the equations be- 
ing solved before we give details on the different solvers 
we have implemented. Section 4 discusses the adap- 
tive mesh refinement strategy and implementation. We 
then move on to describe various test problems for rele- 
vant combinations of solvers, reconstruction schemes in 
one, two and three dimensions, with and without AMR. 
Section 6 presents an application of our code to three- 
dimensional relativistic and supersonic jet propagation 
problem. Section 7 discusses 3d GRB jet simulation. We 
summarize our main conclusions in section 8. 

2. EQUATIONS OF SPECIAL RELATIVISTIC 
HYDRODYNAMICS 

The basic equations of special relativistic hydrodynam- 
ics (SRHD) are conservation of rest mass and energy- 
momentum: 

= 0, (1) 

and 

T% = 0, (2) 

where p is the rest mass density measured in the fluid 
frame, = W(l,v l ) is the fluid four-velocity (assum- 
ing the speed of light c = 1), W is the Lorentz factor, 
v l is the coordinate three-velocity, T^ v is the energy- 
momentum tensor of the fluid and semicolon denotes co- 
variant derivative. 

For a perfect fluid the energy-momentum tensor is 

T^ = phuV+pg^, (3) 

where h = l + e+p/pis the relativistic specific enthalpy, 
e is the specific internal energy, p is the pressure and g^ v 
is the spacetime metric. 

SRHD equations can be written in the form of conser- 
vation laws 



dU dF^ _ 

7)7 + 1^ 7m _ ' 



(4) 



where the conserved variable U is given by 

U = (D,S\S 2 ,S 3 ,t) t , (5) 
and the fluxes are given by 

F j = (Dv j ,S 1 v j ^pS j 1 ,S 2 v j ^pS j 2 ,S 3 v j ^pS j 3 ,S j -Dv j ) T . 

(6) 

lAnile I (|l989l ) has shown that system (|4]) is hyperbolic 
for causal EOS, i.e., those satisfying c s < 1 where the 
local sound speed c s is defined as 



dp 
h [dp 



dp 



(7) 



The eigenvalues and left and right eigenvectors of the 
characteristic matrix dF/dU, whic h are used in some o f 
our numerical schemes, are given bv lDonat etHl (1998). 

The conserved variables U are related to the primitive 
variables by 



D = pW, 
S j =phW 2 v j , 
T = phW 2 -p- pW, 



(8) 
(9) 
(10) 



where j = 1,2,3. The system (j4]) are closed by an equa- 
tion of state (EOS) given by p — p(p,e). For an ideal 
gas, the EOS is, 

p=(r-l)pe, (11) 
where T is the adiabatic index. 

3. NUMERICAL SCHEMES FOR SRHD 

3.1. Time Integration 

We use method of lines to discretize the system (J2J) 
spatially, 



zpx rpx 

r i+l/2J,k r i-l/2,j,k 



dt 



Ax 



F y 

1 i,j + l/2,k 



F y 



Ay 

r i,j,k+l/2 jj,fc-l/2 



Az 



(12) 



L(U), (13) 



where z, j, k refers to the discrete cell index in x, y, z di- 
rections, respectively. F? ±1/2Jtk , i^ i± i /2 , fc and F ij,k±i/2 
are the fluxes at t he cell interface . 

As discussed bv iShu fc OsheTl (fl988h . if using a high 
order scheme to reconstruct flux spatially, one must also 
use the appropriate multi-level total variation diminish- 
ing (TVD) Runge-Kutta schemes to integrate the ODE 
system (p~5|) . Thus we implemented the second and third 
order TVD Runge-Kutta schemes coupled with AMR. 

The second order TVD Runge-Kutta scheme reads, 



uW = U n + A£L(C/ n ), 



(14) 



U n+1 = hj n + hj^ + ^AtL(/7 (1) ), (15) 
and the third order TVD Runge-Kutta scheme reads, 

= U n + AtL(U n ) (16) 
= \u n + + ~AtL(uM) (17) 



U n+l = _ u n 



\uW + \AtL{U^l (18) 



where U n+1 is the final value after advancing one time 
step from U n . 

For an explicit time integration scheme, the time step is 
constrained by the Courant-Friedrichs-Lewy (CFL) con- 
dition. The time step is determined as 



dt = Cmim 



Ax i 



(19) 



where C is a parameter called CFL number and a 1 is 
the local largest speed of propagation of characteristics 
in the direction % whos e explicit expression can be found 
in lDonat et al~l (|1998l ). 
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3.2. Reconstruction Method 

Generally speaking, there are two classes of spatially 
reconstruction schemes (see e.g. LeVeque 2002). One is 
reconstructing the unknown variables at the cell inter- 
faces and then use exact or approximate Riemann solver 
to compute the fluxes. Another is direct flux reconstruc- 
tion, in which we reconstruct the flux directly from fluxes 
at the cell centers. To explore the coupling of different 
schemes with AMR as well as exploring which method 
is most suitable for a specific astrophysical problem, we 
implement several different schemes in both classes. 

To reconstruct unknown variables, we have imple- 
mented piecewise linear method (PLM, Van Leer 1979), 
piecewise parabolic method (PPM, Colella & Woodward 
1984, Marti & Miiller 1996), the third-order convex es- 
sentially non-oscillatory scheme (CENO, Liu & Osher 
1998). These are used to reconstruct the primitive vari- 
ables since reconstructing the conserved variables can 
produce unphysical values in SRHD. Furthermore, un- 
physical values of three-velocities may arise during the 
reconstruction especially for ultrarelativistic flows. So 
we either use v l W to do the reconstruction or we also re- 
construct the Lorentz factor and use it to renormalize the 
reconstructed three-velocity when they are unphysical. 

For direct flux reconstruction, we have implemented 
PLM and the third and fifth order WENO scheme of 
I Jiang fc Shu I (|l996). Direct flux reconstruction us- 
ing WENO was first used t o solve SRHD problems by 
IZhang fc MacFadvenl (|2QQ6j ). They showed that the 
fifth order WENO scheme works well with the third or- 
der Runge-Kutta time integration. In our implemen- 
tation, we followed their description closely. For the 
PLM and CENO sc hemes, we used a generalize d min- 
mod slope limiter (|Kurganov fc Tadmor I I2Q0QT ) . For 
given Vi-i,Vi,v i+ i, v i+1 / 2 = v» + 0.5minmod(6>(^ - 
Vi-i),0.5(vi+i-Vi-i),0(vi+i-Vi)), where 1 < < 2. For 
= 2 it reduces to the m onotonized central-difference 
limiter of I Van Leer I (|l977l ). We found that this gener- 
alized minmod slope limiter behaves much b etter than 
a traditional minmod limiter (|LeVeque Il2002f ) especially 
for strong shear flows. In our calculation, — 1.5 is used 
by default. For t he PPM scheme, we us ed the parame- 
ters proposed by iMarti fc Miiller I (|l996f ) for all the test 
problems. For WENO, we used the pa r amete rs suggested 
in the original paper of iJiang fc Shu I (|l996l ). 

3.3. Riemann Solvers 

In the first class of reconstruction methods, given 
the reconstructed left and right primitive variables 
at interfaces, the flux across each interface is calcu- 
lated by solving the Riemann problem defined by those 
two states. An exact Riemann solver is quite ex- 
pensi ve in SRHD (|Marti fc Miiller I 11994 IPons et al. 1 
2000). Thus we have implemented several approximate 
Riemann solvers including HLL (Harte n et al. I 1 19831 : 
ISchneider et al. I 119931: iKurganov et al. I 120011) . HLLC 
(jToro et al. 1119941 : IMignone fc Bodo I l2QQ5h , local Lax- 
Friedrichs (LLF, Kurga nov fc Tadmor 2 000) and the 
modified Marquina flux (|Alov et al. Ill999h . 

The HLL C scheme is an exte nsion of the HLL solver de- 
veloped by iToro et al71 (|1994f ) for Newtonian flow which 
was extended to t wo-di mensional relativistic flows by 
IMignone fc Bodo I (|2005h . The improvement of HLLC 



over HLL is restoring the full wave structure by con- 
structing the two approximate states inside the Riemann 
fan. The two states can be found by the Rankine- 
Hugoniot conditions between those two states and the 
reconstructed states. With this modification, HLLC in- 
deed behaves better than other Riemann solvers in some 
Id (§ 5.1.7) and 2d (§ 5.2) test problems. But when 
we apply HLLC to three dimensional jet simulation, we 
found that HLLC suffers from the so called "carbuncle" 
artifact we ll-known in th e computational fluid dynamics 
literature (|Quirk I fl994l ). We have used HLLC to run 
many other two-dimensional test problems designed to 
detect the carbuncle artifact and confirmed this short- 
coming. We found that the HLLC solver is unsuitable 
for many multi-dimensional problems. The discussion of 
these problems will be presented elsewhere. In this work, 
we will only apply HLLC to two test problems showing 
that the HLLC solver has less smearing at contact dis- 
c ontinuities than other schem es. 

iLucas- Serrano et al.l (|2004l ) has compared the HLL 
scheme, the LLF scheme and the modified Marquina flux 
formula using Id and 2d test problems. They found those 
three schemes give similar results for all their test prob- 
lems. In our tests, we found similar results. However, 
the modified Marquina flux formula is not as stable as 
HLL in problems with strong transverse flows and LLF is 
more diffusive than HLL. So in the following discussion 
we will only show the results using HLL in most of the 
tests if there is no difference among those three schemes. 

In the following discussion, we will denote a specific 
hydro algorithm by X-Y where X is the flux formula 
and Y is the reconstruction scheme. For example, F- 
WEN05 denotes direct flux reconstruction using fifth 
order WENO. We used the third order Runge-Kutta 
method for all the tests in this work. 

3.4. Converting Conserved Variables to Primitive 
Variables 

Since primitive variables are needed in the reconstruc- 
tion process, after every RK time step, we need to con- 
vert conserved variables to primitive variables. While 
conserved variables can be computed directly from prim- 
itive variables using Eqs. ([8]), ([9]) and ([T0]h the inverse 
operation is not straightforward. One needs to solve a 
quartic equation for the ideal gas EOS and a nonlinear 
equation for more complicated EOS. Iteration methods 
are used even for ideal gas EOS, because co mputing the 
solutio n of a quartic is expensive. Following I Alov et al. I 
(1999), we have used a Newton-Raphson (NR) iteration 
to solve a nonlinear equation for pressure to recover prim- 
itive variables from conserved variables. Typically, the 
NR iteration needs only 2 to 3 steps to converge. 

3.5. Curvilinear coordinates 

We have also implemented cylindrical and spherical 
coordinates following the description of Zhang & Mac- 
Fadyen (2006). This affects three parts of the code. 
Firstly, the geometric factors are incorporated into the 
flux when updating the conserved variables. Secondly, 
there will be geometric source terms. Thirdly, the flux 
correction in AMR (§ 4.4) is modified by geometric fac- 
tors. 

4. ADAPTIVE MESH REFINEMENT 
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4.1. Overview 

Structur ed adaptive mesh refine men t (AMR) was de- 
velope d bv lBerger fc Oligerl (|l984l ) and lBerger fc Colella 
(1989j) to achieve high spatial and temporal resolution 
in regions where fixed grid resolution is insufficient. In 
structured AMR, a subgrid will be created in regions of 
its parent grid needing higher resolution. The hierarchy 
of grids is a tree structure. Each grid is evolved as a sep- 
arate initial boundary value problem, while the whole 
grid hierarchy is evolved recursively. 

renzo is built on top of the AMR framewor k of enzo 
(|Brvan fc Norman 1997a; O'Shealtal] [2004). enzo s 
implementation of AMR follows closely the Berger & 
Colella paper and has been shown to be very efficient for 
very high dynamic range cosmological simulations (see 
e.g. Abel et al. 2002). The pseudocode of the main loop 
for the second order Runge-Kutta method reads, 

EvolveLevel I 

SetBoundaryCondition 
while(fy < ti-i) 

ComputeTimeStep dt\ 
for every grid patch on this level 
Runge — Kutta first step: Eq.(14) 
ComputeFlux 
SweepX, Y, Z 

ChooseHydroAlgorithm 
SaveSubgridFlux 
UpdateConservedVariables 
ConservedToPrimitive 
UpdateTime : t\ = t\ + dt\ 
SetBoundaryCondition(§ 4.2) 
for every grid patch on this level 
Runge — Kutta second step: Eq.(15) 
ComputeFlux 
UpdateConservedVariables 
ConservedToPrimitive 
SetBoundaryCondition 
EvolveLevel I + 1 
UpdateFromFinerGrid (§ 4.4) 
FluxCorrection (§ 4.4) 
RebuildHierarchy for level I 

The RebuildHierarchy function called at the end of ev- 
ery time step is at the heart of AMR. Its pseudocode as 
implemented originally in enzo reads, 

RebuildHierarchy for level I (20) 
for ii eve i = I to MaximumLevel — 1 
for every grid on i level 

FlagCellsForRefinement (§ 4.3) 
CreateSubgrids 

AddLevel (i^vei + 1) 
for every new subgrid 

InterpolateFieldValuesFromParent (§ 4.2) 



CopyFromOldSubgrids 
LoadBalanceGrids 

4.2. Interpolation 

When a new subgrid is created, the initial values on 
that grid are obtained by interpolating spatially from 
its parent grid. In this case, we apply the conservative 
second order order interpolation routine provided by enzo 
to conserved variables. But in this process sometimes the 
interpolated values can violate the constraint (r + D) 2 > 
S 2 + D 2 . If this happens, we will then use first order 
method for that subgrid. 

Before the first Runge-Kutta step for a grid at level I 
(in the following discussion, we use the convention that 
top grid has level 0), we will need the boundary condition 
at time t\ , which is derived by interpolating from its par- 
ent grid. Then at the later steps of Runge-Kutta scheme, 
one needs the boundary condition at time t\ + dt\. Since 
the variables of its parent grid has already been evolved 
to time ti-i + dti-i, which is greater than time t\ + dti, 
we can obtain the boundary conditions at time t\ + dt\ 
for a grid at level I by interpolating both temporally and 
spatially from its parent grid. There are two exceptions 
to this procedure. First, if a cell of fine grid abuts the box 
boundary, then we just use the specified boundary con- 
dition for that cell. Second, if a cell abuts another grid 
at the same level, we copy the value from that grid. Be- 
cause of the above mentioned problem for interpolating 
conserved variables, when interpolating boundary values, 
we apply the second order interpolation to primitive vari- 
ables. Since for ultrarelativistic flows spatially interpo- 
lating three- velocity can lead to unphysical values, we 
also interpolate the Lorentz factor and then use it to 
renormalize the interpolated three- velocity. 

4.3. Refinement Criteria 

In the test problems discussed in section [5j we mainly 
used two general purpose refinem ent criteria that have 
been widely used in AMR code (Zh ang fc MacFadvenl 
l2006h . 

In the first one, we compute the slope 

Si = k+i "7-;' , (2D 

max(|i^|,e) 

where m is typically density, pressure and velocities, e is 
a small number typically taken to be 10 -10 . When Si is 
larger than a minimum slope, typically 1, a cell will be 
flagged for refinement. 
In the second one, for every cell we compute 

E _ \Uj+2 ~ 2Uj + Uj- 2 \ 

\u i+ 2 - Ui\ + \ui - ui- 2 \ + e{\u i+2 \ + 2\ui\ + \ui- 2 \) ' 

(22) 

which is the ratio of the second and first derivatives with 
a safety factor in the denominator. Unless otherwise 
stated, we use e = 0.01. When E{ is larger than a critical 
value E re f, a cell will be flagged for refinement. Typi- 
cally we use E re f =0.8. 

To fully exploit AMR, it is desirable to design more 
specific refinement criterions that are most efficient for a 
specific astrophysical problems. 
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Fig. 1. — Relativist ic blast wave I at t = 0.4 with uniform 
resolution N = 400 for (a) HLL-PLM, (b) HLL-PPM, (c) HLL- 
CENO and (d) F-WEN05. Numerical profiles of density (squares), 
pressure (cross signs) and velocity (diamonds) are shown as well as 
the analytical solution (solid lines). The CFL number used is 0.5. 

4.4. Flux Correction 

When a cell is overlayed by a finer level grid, then the 
coarse grid value is just the conservative average of the 
fine grid values. On the other hand, when a cell abuts 
a fine grid interface but is not itself covered by any fine 
grid, we will do flux correction for that cell, i.e. we will 
use the fine grid flux to replace the coarser grid flux in 
the interface abutting the fine grid, (see Berger & Colella 
1989 for more detailed description of flux correction). For 
this purpose, note that the second order Runge-Kutta 
method can be rewritten as 

U n+1 = U n + ^AtL(U n ) + ^A£L(Z7 (1) ), (23) 

and the third order Runge-Kutta method can be rewrit- 
ten as 

jjn+i = u n Jr I A tL(/7 n ) + -AtL(Z7 (1) ) + -AtL(U i2) ). 
6 6 3 

(24) 

Thus for example, when we do flux correction in the re- 
direction for interface z + 1/2, we will use F^_ 1 ^ 2 (U n )/2-\- 

^+i/ 2 (f /(1) )/2 and F*. 1/2 {U n )/6 + F* +1/2 (U^)/(> + 
2i^ 1 / 2 (k r( ^)/3 to correct the coarser grid conserved 
variables for the second and third order Runge-Kutta 
method, respectively. 

4.5. Parallelism 

renzo uses the enzo parallel framework which uses dy- 
namically load balancing using the Message Passing In- 
terface (MPI) library. At run time, the code will move 
grids among processors according to the current load 
of every processor to achieve a balanced distribution of 
computational load among processors. The computa- 
tional load of a processor is defined as the total number 
of active cells on that processor and level. 

5. CODE TESTS 

5.1. One Dimensional Test 

Relativist ic Riemann problems have analytical solu- 
tions (|Pons et al. I l2QQQh . thus they are ideal for tesing 
SRHD codes. In the following discussion, subscript L 



Fig. 2. — LI errors in rest mass density for the relativistic blast 
wave I problem for six different uniform grid resolutions. The sym- 
bols denote HLL-PLM (plus signs), HLL-PPM (diamonds), HLL- 
CENO (squares) and F-WEN05 (cross signs). The dashed line 
indicates first order of global convergence. 

and R refer to the left and right initial states, respec- 
tively. The initial discontinuity is always at x = 0.5. 
We will report the error between numerical solutions 
and analytical solutions using L\ norm defined as L\ = 
T>i\ui — u(x l )\Ax l where Ui is the numerical solution, 
u (x z ) is the analytical solution and Ax % is the cell width. 

5.1.1. Relativistic Blast Wave I 

This test and the following one are fairly standard and 
all modern SRHD codes can match the analytical solu- 
tion quite well (see Martt & Miiller 2003 for a summary 
of different codes' performance on those two tests). 

The initial left and right states for this problem are 
p L = 13.33, p L = 10.0, vl = 0.0 and p R = lO" 6 , p R = 
1.0, vr = 0.0. The gas is assumed to be ideal with an 
adiabatic index T = 5/3. The initial discontinuity gives 
rise to a transonic rarefaction wave propagating left, a 
shock wave propagating right and a contact discontinuity 
in between. This problem is only mildly relativistic with 
a post-shock velocity 0.72 and shock velocity 0.83. The 
results using four hydro solvers are shown in Fig. [TJ The 
CFL number used is 0.5. The L\ errors are shown in 
Fig. [2j from which we can see the four schemes behave 
essentially identical for this problem. The order of global 
convergence rate is about 1 for all four schemes, which is 
consistent with the fact that there are discontinuities in 
the problem. 

5.1.2. Relativistic Blast Wave II 

The initial left and right states for this problem are 
p L = 1000.0, pL = 1.0, vl = 0.0 and p R = 10~ 2 , p R = 
1.0, vr = 0.0. The gas is assumed to be ideal with an 
adiabatic index T = 5/3. This test is more relativistic 
than the previous one. While the wave structure is the 
same, the thermodynamically relativistic initial left state 
gives rise to a relativistic shock propagating at a Lorentz 
factor W ~ 6 and a very thin dense shell behind the 
shock with width w 0.01056 at t = 0.4. The CFL number 
used is 0.5. The results using four hydro algorithms are 
shown in Fig. [3j The L\ errors are shown in Fig. 0J 
We can see that for this problem PPM and WENO have 
smaller L\ errors than PLM and CENO. This is due to 
their ability to better resolve the thin shell. 
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Fig. 3. — Relativist ic blast wave II at t = 0.4 with uniform 
resolution N = 400 for (a) HLL-PLM, (b) HLL-PPM, (c) HLL- 
CENO and (d) F-WEN05. Numerical profiles of density (squares), 
pressure (cross signs) and velocity (diamonds) are shown as well as 
the analytical solution (solid lines). The CFL number used is 0.5. 



Fig. 6. — LI errors in rest mass density for the planar jet problem 
for six different uniform grid resolutions. The symbols denote HLL- 
PLM (plus signs), HLL-PPM (diamonds), HLL-CENO (squares) 
and F-WEN05 (cross signs). The dashed line indicates first order 
of global convergence. 





Fig. 4. — LI errors in rest mass density for the relativistic blast 
wave II problem for six different uniform grid resolutions. The 
symbols denote HLL-PLM (plus signs), HLL-PPM (diamonds), 
HLL-CENO (squares) and F-WEN05 (cross signs). The dashed 
line indicates first order of global convergence. 
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Fig. 5. — Planar jet problem at t = 0.4 with uniform resolution 
N = 400 for (a) HLL-PLM, (b) HLL-PPM, (c) HLL-CENO and (d) 
F-WEN05. Numerical profiles of density (squares), pressure (cross 
signs) and velocity (diamonds) are shown as well as the analytical 
solution (solid lines). The CFL number used is 0.5. 



Fig. 7. — Blast wave with transverse velocity problem I at t = 0.4 
with uniform resolution N = 400 for (a) HLL-PLM, (b) HLL-PPM, 
(c) HLL-CENO and (d) F-WEN05. Numerical profiles of density 
(squares), pressure (cross signs) and velocity (diamonds) are shown 
as well as the analytical solution (solid lines). The CFL number 
used is 0.5. 

5.1.3. Planar Jet Propagation 

The initial left and right states for this test are Pl = 
1.0, Pl = 1.0, v L = 0.9 and p R = 10.0, p R = 1.0, 
vr = 0.0. The gas is assumed to be ideal with an adia- 
batic index T = 4/3. This test mimics the interaction of 
a planar jet head with the ambient medium. The decay 
of the initial discontinuity gives rise to a strong reverse 
shock propagating to the left, a forward shock propagat- 
ing to the right and a contact discontinuity in between. 
The results are shown in Fig. [5l The CFL number is 
0.5. The L\ errors are shown in Fig. [6l As can be 
seen in Figs. [5] & El for this problem PPM and WENO 
behave better than PLM and CENO: there is almost no 
oscillation behind the reverse shock and they capture the 
contact discontinuity with fewer cells. Especially PPM 
captures the contact discontinuity with only 4 cells and 
has the smallest L\ error. 

5.1.4. Blast Wave with Transverse Velocity I 
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Fig. 8. — LI errors in rest mass density for the blast wave with 
transverse velocity problem I for six different uniform grid reso- 
lutions. The symbols denote HLL-PLM (plus signs), HLL-PPM 
(diamonds), HLL-CENO (squares) and F-WEN05 (cross signs). 
The dashed line indicates first order of global convergence. 

For this and the following two problems, we will con- 
sider non-zero transverse velocities in the initial states. 
The initial state is identical to blast wave problem II ex- 
cept the presence of transverse velocities. Those prob- 
lems were first discussed analytically by iPons et al. I 
(|2000f ) . Since then various groups have shown that when 
transverse velocities are non-zero, in some cases those 
problems become very difficult to solve numerically un- 
less very high spatial resolution is used ([M ignone et al.l 
120051 : IZhang fc MacFadvedl2006t iMorsonv et al. Il2006f ). 
In realistic astrophysical phenomena transverse veloci- 
ties are usually very important (see e.g. Aloy & Rezzolla 
2006), thus solving those problems accurately is of great 
importance. 

As an easy first case, we will consider non-zero trans- 
verse velocity only in the low pressure region. The ini- 
tial left and right states are pl = 1000.0, pL = 1-0, v x l = 
0.0, v yL 0.0 and p R 10~ 2 , p R 1.0, v xR 0.0, v yR 
0.99. The gas is assumed to be ideal with an adiabatic 
index T = 5/3. The results are shown in Fig. [7| The 
CFL number is 0.5. The L\ errors are shown in Fig. 
[H We can see that all four hydro algorithms behaves 
similarly well, except that PLM and CENO shows some 
small oscillations around the contact discontinuity. 

5.1.5. Blast Wave with Transverse Velocity II 

Next, we consider non-zero transverse velocity in the 
high pressure region. In this case, the problem becomes 
more difficult to solve num erically (|Mignone et al.l 2005; 
IZhang fc MacFadvenll2006h . The initial left and right 
states for this problem are Pl — 1000.0, pL = 1.0, v x l = 
0.0, v yL 0.9 and p R 10~ 2 , p R 1.0, v xR 0.0, v yR 
0.0. The gas is assumed to be ideal with an adiabatic 
index T = 5/3. The high pressure region is connected to 
the intermediate state by a rarefaction wave. Since the 
initial normal velocity in the high pressure region is zero, 
the slope of the adiabat increases rapidly with transverse 
velocity, thus a large initial transverse velocity will lead 
to a small intermediate pressure and a small mass flux. 

The results using a uniform grid and two AMR runs 
are shown in Fig. [9l The hydro solver used for this fig- 
ure is HLL-PLM. The CFL number is 0.4. The L x error 
are shown in Fig. [TUl It can be seen that for the run 
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Fig. 9. — Blast wave with transverse velocity problem II at t = 
0.4 with uniform resolution N = 400 (top), 4 refinement levels with 
a refinement factor of 2 (middle, equivalent resolution 6400) and 3 
refinement levels with a refinement factor of 4 (bottom, equivalent 
resolution 25600). The hydro solver used is HLL-PLM. Thick lines 
are numerical solution while thin lines are analytical solutions. The 
CFL number used is 0.4 




Fig. 10. — LI errors in rest mass density for the blast wave with 
transverse velocity problem II for seven different equivalent grid 
resolutions using AMR. The solid lines with square signs are for 
HLL-PLM. The dashed line indicates first order of global conver- 
gence. 

with 400 uniform grid cells, the n umerical solution is in - 
adequate, as previously found by iMignone et all ([2005). 
This is mainly due to the poor capture of the contact 
discontinuity. We have tried to run this problem with 
various algorithms but only obtained accurate solutions 
by dramatically increasing the resolution. 

5.1.6. Blast Wave with Transverse Velocity III 

Now we introduce transverse velocity in both region. 
The initial left and right states for this problem are 
p L 1000.0, pL = 1.0, v x l 0.0, Vy L 0.9 and 
Pr = lO~ 2 ,p R = l.0,v xR = 0.0, v yR = 0.9. The gas is 
assumed to be ideal with an adiabatic index T = 5/3. 
This problem is more difficult than the previous one 
due to the formation of an extremely thin shell between 
the rarefaction wave tail a nd the contact discontinuity 
(jZhang fc MacFadvenll2006h . 
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Fig. 11. — Blast wave with transverse velocity problem III at t = 
0.4 with uniform resolution N = 400 (top), 4 refinement levels with 
a refinement factor of 2 (middle, equivalent resolution 6400) and 4 
refinement levels with a refinement factor of 3 (bottom, equivalent 
resolution 25600). The hydro solver used is HLL-PLM. Thick lines 
are numerical solutions while thin lines are analytical solutions. 
The CFL number used is 0.5 



TABLE 1 
Equivalent resolution and the 
actual number of cells used in blast 

Wave problems with Transverse 
Velocity II (BT II) and III (BT III). 



Equivalent Resolution 


BT II 


BT III 


400 


400 


400 


800 


448 


421 


1600 


455 


442 


3200 


470 


468 


6400 


501 


473 


12800 


518 


476 


25600 


562 


520 




Fig. 12. — LI errors in rest mass density for the blast wave with 
transverse velocity problem III for seven different equivalent grid 
resolutions using AMR. The solid lines with square signs are HLL- 
PLM. The dashed line indicates first order of global convergence. 




The results with a uniform grid and two AMR runs are 
shown in Fig. [TH The hydro solver used for this run is 
HLL-PLM. The CFL number used is 0.5. The L\ errors 
are shown in Fig. fT2l 

Table 1 shows the equivalent resolution and the ac- 
tual number of cells used for this and the previous tests. 
It can be seen for the highest resolution calculation our 
code uses about four hundred times less grid cells than 
the corresponding uniform grid calculation. Thus AMR 
allows us to achieve very high resolution while signifi- 
cantly reducing the computational cost. 

5.1.7. Jet-Cocoon interaction 

For this test we set up a one-dimensional Riemann 
problem that mimics the interaction of jet with an over- 
pressured cocoon. The initial left and right states for 
this problem are pl = 0.00017, pL = 0.01, v x l = 0, v v l = 
0.99 and p R = 0.017, p fl = 0.1, = 0,^ = 0. The 
gas is assumed to be ideal with T = 5/3. Those val- 
ues mimic th e conditions of the jet-cocoon boundary in 
model C2 of iMarti etafl (|l997f ). The result using four 
different Riemann solver with PLM on uniform grid are 
shown in Fig. [T3l It can be seen that solutions use 
HLL, LLF and direct flux reconstruction have large pos- 
itive fluctuations in the normal velocity at the rarefaction 



Fig. 13. — Jet-cocoon interaction problem at t = 0.4 for (from 
top to bottom): HLL-PLM, LLF-PLM, HLLC-PLM and F-PLM. 
Thick lines are numerical solution while thin lines are analytical 
solutions. The grid resolution is uniform TV = 400. The CFL 
number is 0.5. 

wave. It is interesting to note that only HLLC does not 
suffer from this shortcoming which is probably due to the 
ability of HLLC to resolve the contact discontinuity com- 
pared to the other Riemann solvers in the code. If those 
fluctuations also happen in higher dimensional jet simu- 
lation, then one would expect that the normal velocity 
fluctuation seen in this test would lead to an artificially 
extended cocoon. 

In Fig. M the result of using HLL-PLM and HLLC- 
PLM with AMR is shown. It can been seen that the fluc- 
tuation in the HLL scheme becomes smaller with higher 
resolution. Fig. [15] shows the L\ error for those two 
schemes with different levels of refinement. 

5.2. Two Dimensional Test: Shock Tube 

To test the code in higher dimension, we first study 
the two-dimensional shoc k tube problem suggested by 
iDel Zanna fc Bucciantini I (|2002l ) and latter also used 
by various groups (see e.g. Zhang & MacFadyen 2006, 
Mignone & Bodo 2005a, Lucas-Serrano et al. 2004). This 
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Fig. 14. — Jet-cocoon interaction problem at t = 0.4 for HLL- 
PLM and HLLC-PLM. The top two are HLL-PLM and HLLC- 
PLM with 4 refinement levels and a refinement factor of 2 (equiva- 
lent resolution 6400). The bottom two are HLL-PLM and HLLC- 
PLM with 4 refinement levels and a refinement factor of 4 (equiv- 
alent resolution 102400). Thick lines are numerical solutions while 
thin lines are analytical solutions. The CFL number is 0.5. 




Fig. 15. — L\ error for the jet-cocoon interaction problem for 
eight different grid resolutions using AMR. The square signs are for 
HLL-PLM while the diamonds are HLLC-PLM. The dashed line 
indicates first order of global convergence. 

test is done in a two-dimensional Cartesian box divided 
into four equal- area-constant states: 

(p,v x ,v y ,p) NE = (0.1,0,0,0.01), (25) 

(p,v x ,v y ,p) NW = (0.1,0.99,0,1), 

(p,v x ,v y ,p) sw = (0.5,0,0,1), 

(p,v x ,v y ,p) SE = (0.1,0,0.99,1), 

where NE means northeast corner and so on. The grid is 
uniform 400 x 400. The gas is assumed to be ideal with 
an adiabatic index T = 5/3. We use outflow boundary 
conditions in all four directions and the CFL number is 
0.5. 

The results are shown in Fig. [16] for four schemes. 
This problem does not have analytical solutions to com- 
pare with, but comparing our result with other groups' 
result shows good agreement. The cross in the lower left 
corners of (a) HLL-PLM and (d) F-WENO is a numer- 




0,0 0,2 0,4 0,6 o,; 



1,0 0,2 0,4 0,6 0,3 1,0 



Fig. 16. Shock tube problem at t = 0.4 for (a) HLL-PLM, (b) 
HLL-PPM, (c) HLLC-PLM and (d) F-WEN05. Thirty equally 
spaced contours of the logarithm of density are plotted. The CFL 
number is 0.3. 
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Fig. 17. — One-dimensional cut of the relativistic spherical shock 
reflection problem at t = 0.4 for density (top), pressure (middle) 
and radial velocity (bottom). The solid line is the analytical solu- 
tion. Square signs, plus signs and dashed lines are cut though x 
axis, y axis and diagonal direction, respectively. The resolution is 
a uniform grid with 128 3 zones. LLF-PLM are used with a CFL 
number of 0.1. It can be seen that the fluctuations in density and 
pressure at the sphere center present in many previous codes are 
absent in our results. 

ical artifact due to the inability to maintain a contact 
discontinuity perfectly, which are absent in the results 
us ing the HLLC solv e r (c). This agreed with the result 
of iMignone fc Bodo I (|2005f ) that the HLLC solver be- 
haves better in this problem than other Riemann solvers 
because of its ability of resolve contact discontinuities. 

5.3. Three Dimensional Tests 
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5.3.1. Relativistic Spherical Shock Reflection 

We first study a test problem with uniform grid 
in Cartesian coordinate that has been used by sev- 
eral groups to test the symme tric propert ies of 

a three-dimensiona l SRHD code (|Aloy et al. I 1 19991 : 
IMignone et al.ll2005l ). The initial setup of this test con- 
sists of a cold spherical inflow with initially constant den- 
sity po and constant velocity v\ colliding at the box cen- 
ter. This problem is run using three-dimensional Carte- 
sian coordinates so it allo ws one to evaluate the symme- 
try p r operties of the code ([Marti et al.ll 19971 : lAlov et al~l 



1999; IMignone et all |2005[ ) . When gas collides at the 
center, a reflection shock forms. Behind the shock, the 
kinetic energy will be converted completely into internal 
energy. Thus the downstream velocity V2 is zero and the 
specific internal energy is given by the upstream specific 
kinetic energy 

e 2 = Wi - 1. (26) 

Using the shock jump condition, the compression ratio 
p2 1 Pi and the shock velocity can be found to be 



P2 

Pi 



r- i r- i 
_ (r-i)Wih 



^2 



1 



(27) 
(28) 



In the unshocked region (r > V 3 t), the gas flow will 
develop a self-similar density distribution, 



Pi 



1 + ** 
r 



po- 



rn 



The initial state are pi = 7 x 10 -6 , po = 1, v\ = —0.9. 
We chose a small value for pressure because a grid-based 
code cannot handle zero pressure. A CFL number 0.1 
is us ed for this problem, as other groups (|Alov et al. I 
11999( 1. We chose to use LLF-PLM for this problem be- 
cause this turns out to be the most stable solver for this 
problem. Fig. [T7l shows the one-dimensional cut though 
axis and diagonal direction and Fig. [18] shows a contour 
through z = 0.5 plane, both at t = 0.4. It can be seen 
from those plots that our code keeps the original spher- 
ical symmetry quite well. Since in a Cartesian box the 
simple outflow boundary condition is inconsistent with 
the initial spherical inflow setup, we evolve this problem 
only to t = 0.4, at which point all the mass in the origi- 
nal box has just entered the shocked region (see Fig. fTTj) . 
After that time, the evolution would be affected by the 
unphysical boundary condition. 

5.3.2. Relativistic Blast Wave I 

In this test, we study a spherical blast wave in three- 
dimensional Cartesian coordinates. There is no analyti- 
cal solution for this problem. Thus for the sake of com- 
parison, we set up the same problem as other groups 
(iDel Zanna fc BucciantnTil 120021 : IZhang fc MacFadvenl 
2006). The center of the blast wave source is located 
at the corner (0,0,0) of the box. The initial conditions 
are 



(p,v r ,p) = (1,0, 1000) if r<0.4, 
(p,v r ,p) = (1,0,1) if r>0.4, 
where r is the distance to the center (0, 0, 0). 



(30) 




0.0 0.2 0.4 0.6 0.8 1.0 

Fig. 18. — Density contour through z = 0.5 of of the relativistic 
spherical shock reflection problem at t = 0.4. It can be seen that 
the symmetry of the initial condition is preserved rather well up to 
this time. After this time, the unphysical boundary condition will 
begin to affect the subsequent evolution. 




Fig. 19. — One-dimensional cut of the three-dimensional blast 
wave problem I at t = 0.4 for density (top), pressure (middle) and 
velocity (bottom). The triangle signs and cross signs are cuts along 
the x axis and diagonal direction, respectively. The solid line is a 
one-dimensional simulation run using spherical coordinates with 
4000 uniform zones. F-PLM is used for this calculation. A top 
grid of 128 3 zones with two levels of refinement and a refinement 
factor of 2 is used. The CFL number is 0.1. 

An ideal gas with an adiabatic index of T — 5/3 is 
assumed. The left boundaries at x, y, z directions are 
reflecting while others are outflow. We use a top grid 
of 128 3 zones with two more levels of refinement and 
a refinement factor of 2 for this calculation (equivalent 
resolution 512 3 ). F-PLM is used for the result shown and 
the CFL number is 0.1. 

The results are given in Fig. fT9l which shows the cut 
along x-axis and diagonal direction. For comparison, 
we run a high resolution one-dimensional simulation us- 
ing spherical coordinates. The three-dimensional run in 
Cartesian coordinates agrees with the one-dimensional 
high resolution run. Furthermore, it can be seen that the 
spherical symmetry of the initial condition is preserved 
rather well in the three-dimensional Cartesian case. 

5.3.3. Relativistic Blast Wave II 

Finally we study another blast wave problem for which 
the center of the blast wave source is located at the box 
center. This problem also does not have analytical solu- 
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Fig. 20. — Three-dimensional blast wave problem II at t = 0.12. 
All points in a slice at y = 0.5 are plotted as a function of the 
distance to the center (0.5,0.5,0.5). The upper panel shows the 
laboratory density D and the lower panel shows the Lorentz factor 
W. HLL-PLM is used for this calculation. A top grid of 64 3 zones 
with four levels of refinement and a refinement factor of 2 is used. 
The CFL number is 0.5. 

tion but it has been studied by ( Hug hes" et al.ll2QQ2l ) so 
our result can be compared to theirs. The initial condi- 
tions are 

(p,v r ,p) = (1,0, 10 4 ) if r<0.05, (31) 
(p,v r ,p) = (0.1,0, 10) if r>0.05. 

An ideal gas EOS with an adiabatic index of T = 4/3 is 
used. We s top the run at t = 0. 12, roughly the same end- 
ing time of lHughes et al.l (|2002f ). A top grid of 64 3 zones 
with four levels of refinement and a refinement factor 
of two is used (equivalent resoluti on 1024 3 ). Thus ou r 
resolution is roughly 1.5 times of iHughes et al.l ([2002) . 
We used HLL-PLM and a CFL number of 0.5 for this 
calculation. 

Fig. [20] plots the numerical solution for all cells cen- 
tered on the highest level in the two dimensional slice at 
y = 0.5 at t = 0.12 as a function of radius from the cen- 
ter (0.5,0.5,0.5). The position and amplitu de of the high 
densit y shell agrees with the calculation of lHughes et al.l 
(2002). And it can be seen that the spherical symmetry 
is preserved rather well in our code. 

6. ASTROPHYSICAL APPLICATION I RELATIVISTIC 
SUPERSONIC JET PROPAGATION 

Having validated our code using various test problems, 
we can apply it to astrophysical problems. We will study 
two typical astrophysical relativistic flow problems in this 
work, AGN jets and GRB jets. Both topics have been 
studied extensively with 2d simulations before, but very 
few 3d calculations have been done in both cases. Con- 
sequently we will focus on 3d simulations here. 

In this section, we study a relativistic supersonic jet 
in three dimensions. We set up the problem us i ng; th e 
same parameters as model C2 of ([Marti et al.l [l997). 
This model has also b een s tudied in two dimensions by 
Zhan g fc MacFadyenl (|2006l ) and in three dimensions by 
Alov et al. I (|l999h . The jet parameters are p& = 0.01, 
Vb = 0.99 and = 0.02. The jet has a classical 
Mach number M5 = Vb/c s = 6 so the pressure is 
Pb = 0.000170305. The parameters for the medium are 
p rn = 1.0, v m = and p rn = p^. The EOS is assumed 
to be ideal with T = 5/3. The jet is injected from the 
low-z boundary centered at (0.5,0.5,0) with radius ?v 



Outflow boundary conditions are used at other part of 
the boundary. 

Figure [21] shows the result at t — QOr^/c for the 
three runs: HLL-PPM with three levels of refinement 
(HLL-PPM-13), HLL-PLM with three levels of refine- 
ment (HLL-PLM-13) and HLL-PLM with four levels of 
refinement (HLL-PLM- 14). The top grid resolution is 
64 3 zones. Thus the first two runs have an equivalent 
resolution of 512 3 zones while the last one has 1024 3 
zones. For the first two, turbulence in cocoon is not 
fully developed so the cocoon is still symmetric even in 
3d. The HLL-PPM-13 run has slightly more turbulent 
cocoon due to the higher spatial reconstruction order 
of PPM. Thus the HLL-PPM-13 jet propagates slightly 
slower than the HLL-PLM-13 jet. On the other hand, for 
the HLL-PLM-14 jet, the resolution is 20 cells per beam 
radius, comparable to the resolut i on us ed in the two- 
dimensional study by iMarti et al.l (|l997l ). The cocoon 
turbulence is much more dev eloped in this case, as in 
the 2d case ([Marti et al.l fl997). Consequently, the HLL- 
PLM-14 jet propagates slower than the two lower resolu- 
tion ones. Furthermore, the HLL-PLM-14 case does not 
show axisymmetry because instability quickly develops 
in the lateral motion and consequently lateral motion 
also becomes turbulent. 

Since we found the jet-cocoon structure differs signifi- 
cantly at higher resolution run and consequently the jet 
propagation speed decreases, we conclude that even at 
1024 3 effective resolution of our three dimensional jet 
simulations the correct solution remains elusive. More- 
over, different solvers give disparate answers. 

7. ASTROPHYSICAL APPLICATION II. THREE 
DIMENSIONAL SIMULATION OF COLLAPSAR JETS 

The idea that "long-soft" gamma-ray bursts are associ- 
ated with the deaths of massive stars has been supported 
by various observations recently (see e.g. Woosley & 
Bloom (2006) for a recent review). So it is of great inter- 
est to calculate how a relativistic jet propagate through 
and break out of a massive sta r. There has been exten- 
sive calculations in 2d recently (IZhang et al.ll2003l 120041 : 
iMizuta et aI1l2006l : iMorsonv et al. 1120061 ). But very few 
3d calculation have been reported in the literature so far. 
Since there are significant turbulent motion and mixing 
in the jet-cocoon system in 2d calculation, it is important 
to model this process in 3d. 

We use the projenitor HE16A of Wo oslev fc H eger 
(2006), which is a stripped-down helium core with ini- 
tial mass 16 M . Fig. 1221 shows the density and pressure 
profiles for the progenitor model. The radius of the star 
at onset of collapse is 3.86 x 10 10 cm. In our simulation, 
the mass inside 5 x 10 9 cm is removed and replaced by a 
point mass of 9.3 M . The jet has power E = 3 x 10 50 
erg s _1 , initial Lorentz factor To = 5, the ratio of its 
total energy (excluding rest mass energy) to its kinetic 
energy fo = 40. This corresponds to a a jet with initial 
density 5.937 x 10 -3 g cm -3 and pressure 4.156 x 10 19 
erg cm -3 . The jet is injected parallel to the z axis with 
an initial radius 8.73 x 10 8 cm. Those m odel parame- 
ters a re similar to p revious 2d calculations ([Zhang et al.l 
1200a IMizuta et al.ll2006h . 

We use a simulation box of 3.2 x 10 11 cm in order to 
follow the propagation of the jet after break out. We use 
HLL-PLM as this should be the most stable and reliable 
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Fig. 21. — Density slice through y = 0.5 for three dimensional relativistic jet at t=60rb/c for HLL-PPM with three levels of refinement 
(left), HLL-PLM with three levels of refinement (middle) and HLL-PLM with four levels of refinement (right). The top grid resolution is 
64 3 zones and a refinement factor of 2 is used. Thus the left and middle panels have an equivalent resolution of 512 3 zones while the right 
one has 1024 3 zones. The CFL number is 0.4. 

percent of jet material, we flag that cell for refinement. 
This ensures that we also have high resolution when mix- 
ing between jet material and star material happens. 

The results for those two runs are show in Fig. [23j It 
can been seen that the high resolution run gives qualita- 
tively different jet dynamics. While in the low resolution 
run the jet break out of the star successfully, in the high 
resolution run the jet head bifurcate at the stellar edge. 
Interestin gly, this behavior has also been seen in 2d cal- 
culation (Mo rsonv et al. II2QQ6I ). But since we did not 
get convergent behavior so far, we cannot conclude at 
this stage whether this behavior is physical or purely nu- 
merical. However, it is safe to conclude that much higher 
resolution will be needed to model the jet break out in 
3d. 

8. CONCLUSIONS AND DISCUSSIONS 

In this paper, we have described a new code that 
solves the special relativistic hydrodynamics equations 
with both spatially and temporally adaptive mesh refine- 
ment. It includes direct flux reconstruction and four ap- 
proximate Riemann solvers including HLLC, HLL, LLF 
and modified Marquina flux formula. It contains sev- 
eral reconstruction routines: PLM, PPM, third order 
CENO, and third and fifth order WENO schemes. A 
modular code structure makes it easy to include more 
physics modules and new algorithms. 

From our test problems and two astrophysical applica- 




FlG. 22. — Density profile and pressure profile of the HE16A 
progenitor model adopted in our GRB jet calculation. 



scheme to carry out resolution study. The top grid reso- 
lution is 128 3 . We run the simulation using both 4 and 5 
refinement levels, which correspond to resolution of 5.6 
and 11 cells per jet beam radius, respectively. In order 
to always have high resolution for the jet material, we 
designed a color field refinement strategy in addition to 
the standard refinement criterion designed for disconti- 
nuities. More specifically, we use two color fields to keep 
track of the injected jet material and the stellar material. 
Those two color fields give us a fraction of jet material at 
every cell. Then whenever a cell contains more than 0.1 
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Fig. 23. — Density slice through y = 0.5 for three dimensional relativistic GRB jet at t=10.99 s for HLL-PLM with four levels of 
refinement (left) and five levels of refinement (right). The top grid resolution is 128 3 zones and a refinement factor of 2 is used. Thus the 
left and middle panels have an equivalent resolution of 2048 3 zones while the right one has 4096 3 zones. The CFL number is 0.4. 



tions, it is clear that relativistic flow problems are more 
difficult than the Newtonian case. One key reason is that 
in the presence of utrarelativistic speed, nonlinear struc- 
tures such as shocked shells are typically much thinner 
and thus requires the use of very high spatial resolu- 
tion. SRHD problems also become difficult to solve ac- 
curately when significant transverse velocities are present 
in the problem as we have shown using several one dimen- 
sional problems. One reason for this difficulty is that in 
SRHD velocity components are coupled nonlinearly via 
the Lorentz factor. In studying astrophysical jet prob- 
lems, we have demonstrated the need of both high reso- 
lution achievable only through AMR and careful choice 
of hydrodynamic algorithms. In addition to validate our 
AMR code, the most important implications of the cal- 
culations we have done is that in relativisitic flow simu- 
lations, resolution studies are crucial. 
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